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The density of states for the three-dimensional Ising model is calcu- 
lated with high-precision from multicanonical simulations. This allows 
us to estimate the leading partition function zeros for lattice sizes up 
to L = 32. Combining previous statistics for smaller lattice sizes, we 
have evaluated the correction to scaling and the critical exponent v 
through out an analysis of a multi-parameter fit and of Bulirsch-Stoer 
(BST) extrapolation algorithm. The performance of the BST algo- 
rithm is also explored in case of the 2D Ising model, where the exact 
partition function zeros are known. 
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1 Introduction 



In recent years there has been a persistent interest in obtaining accu- 
rate estimates of critical parameters of the three-dimensional (3D) Ising-like 
systems through high performance simulations and perturbative expansions 
0, |, |, g, |, |, 0, |, |, 0. Our aim here is to enlarge the knowledge about 
critical behavior of Ising-like 3D systems by calculating not only the critical 
exponent of the correlation length u, but also by tackling the much harder 
problem of calculating the first correction to scaling w. 

A common way to extract information on phase transitions from Monte 
Carlo simulations is by means of finite-size scaling (FSS), for instance by ana- 
lyzing the partition function zeros in the complex temperature plane |Lj] , [12| . 
This approach is not restricted to Ising-like systems and was recently even 
used to study structural transitions in bio- molecules [DJ. However, separa- 
tion of universality classes can be tricky for 3D systems, and requires high- 
precision data. An important tool for obtaining such data was provided by 



Ferrenberg and Swendsen [TJj who revived reweighting techniques introduced 
by Salsburg et al. [[y| forty years ago. In the sequence, multiple-histogram 



[|T6| , |17|j and multicanonical simulations |Tj| were proposed for a reliable nu- 
merical determination of the density of states and extensively checked in 
two- dimensions where exact results are available. Only by combining one 
of these sophisticated new simulation techniques (exhaustive multicanonical 
Monte Carlo simulations for L < 32) with Itzykson's finite-size scaling rela- 
tion for the partition function zeros we were able to obtain the results 
presented in this article. 

Before proceeding further, we give the outline of the paper. In the next 
section we describe the numerical evaluation of the partition function. Sec- 
tion 3 is concerned with the crucial task of data analysis. There we calculate 
the critical exponent v and the correction to scaling. The exponent v was 
first obtained by a four-parameter fitting. Next, we used the Bulirsch-Stoer 



algorithm [19], g0| to extrapolate results from finite lattices to the thermo- 
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dynamical limit. That algorithm has a free-parameter w which is exactly 
the wanted correction to scaling. We further corroborate on results recently 
obtained by other techniques. Finally, we test the performance of our ap- 
proach, and specially the BST algorithm, for the 2D Ising model where the 
exact zeros are known. 



2 Multicanonical simulation and partition func- 
tion zeros 



In the multicanonical algorithm JTjl conformations with energy E are as- 
signed a weight w mu (E) oc l/n(E). Here, n(E) is the density of states. A 
simulation with this weight will lead to a uniform distribution of energy: 

P mu {E) oc p(E) w mu (E) = const . (1) 

This is because the simulation generates a ID random walk in the energy, 
allowing itself to escape from any local minimum. Since a large range of 
energies are sampled, one can use the reweighting techniques jllj to calculate 
thermodynamic quantities over a wide range of temperatures by 

/ dx A(x) w~\E(x)) e~ mx) 

<A> T = J — r , (2) 

dx w-\E{x)) e~ mx) 



where x stands for configurations. 

It follows from equation [I] that the multicanonical algorithm allows us to 
calculate estimates for the spectral density: 

p(E) = P rnu (E)w-l(E) . (3) 

We can therefore construct the partition function of the three-dimensional 
Ising model 

Z((3)=Y,P(E)u E , (4) 

E 
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where we define u = e -4 ^. 

In the present article, the estimates of the partition function rely on 
averages over NRUN simulations of NSWEEP Monte Carlo updates for 3D 
Ising models of linear size L. In order to allow the system to thermalize, 
additional sweeps in the beginning were performed and discarded. Table 1 
lists the respective values for NSWEEP and NRUN. 

Once, we have calculated reliable estimates for the partition function we 
can calculate the zeros. In the polynomial form, equation (|4]) has a large 
number of coefficients for large lattice sizes. For this reason we apply the 
method presented in |2l[] to obtain the leading complex zeros u\{L). In table 



1 we collect the so calculated zeros as obtained from our multicanonical 
simulations. 

3 Finite size scaling analysis 

The standard FSS approach for the zero u\(L) closest to the real positive 
axis neglects, for sufficiently large L, corrections to scaling [ |Tl"|j , 

u1(L) = u c + AL- 1/v [l + 0{L- W )] . (5) 

Since we are interested in analysing v in function of the correction to scaling 
exponent w, we follow ref. [17], ^| and fit our data to a four parameter 
scaling relation, 

\u\{L) - u c \ = ai L~ 1/u + a 2 L- a3 . (6) 

However, due to the small number of available lattices, such a multi-parameter 
fit is not appropriate. To circumvent this difficulty we proceed including 
statistics already available for smaller lattice sizes. In ref. |p.7| , table 1, zeros 
for 3D Ising model are presented for L — 6,8, 10 and 14, together with data 



from ref. |22| for L = 3,4,5,6,8 and 10. In addition, we decided to take 
into account that information and our results in table 1, weighted with the 
corresponding available precision. The final estimates for each lattice is a 
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linear combination of available data with normalized weight factors which 
are taken as the reciprocals of the corresponding empirical variances for each 
data ^3 . 



In figure 1 we show the corresponding fit for all lattice sizes, whose pa- 
rameters are found by monitoring the goodness of fit Q . Here we choose 
the recently obtained critical value u c = 0.41204684(25) |§, to apply the 
least-square method to equation @, although the fit is not very sensible 
to the precision in the value of critical temperature u c . The error bar of 
those data is included, but it is hardly seen in that scale. The goodness 
of fit (Q = 0.89) reveals a very good agreement with the data. We obtain 
v = 0.62853(35) and a 3 = 4.861(84). 

If we discard the smallest sizes L = 3,4 and 5, we obtain v = 0.6280(15) 
(Q = 0.84), corresponding to y t = \ jv = 1.5924(38), which is in remarkable 
agreement with previous results [fj, [7|, |8|. However, this fit is less stable with 
relation to the parameters in the second term. This is because of the presence 
of rather large lattice sizes. 

3.1 Correction to scaling from RG transformation 

In this section we now evaluate the correction to scaling by a method which 
is similar to the "finite size phenomenological renormalization group (RG)" 
analysis by Binder (which in turn was based on Nightingale's finite size 
RG transformation for the correlation length Our method to evalu- 



ate w was previously used J21| to analyse the 2D Ising model, and it is briefly 
recalled here. 

In order to consider the scaling relation for the longitudinal correlation 
length £l(/3), we assume that the system is of finite length scale L in one 
direction and infinite in all other directions. The standard expression for the 
correlation exponent v is now given by |27], |28| 



1 , (dtu/dp\ n ,L\ 
H = In . ' /ln(— ). (7) 



5 



This expression is obtained from a linearization around the fixed critical 
point (3 C for fixed scaling transformation L — > L'. The scaling equation for 
the finite size longitudinal correlation length is given by 

Cl = L Yg( (P — (3 c )L l l v , hL VH ,uL ys ) . (8) 

This differentiable equation includes corrections due to the leading bulk irrel- 
evant scaling field u with exponent y$ < 0, and a magnetic field dependence 
for the sake of completeness. 

From equation (|7|) and equation (H) one obtains, for h = 0, 

1 1 L' ys - L y3 , L' 2y3 - L 2y3 

a . 777777 + &o ; fTllT . + ... (9) 



v \n(L'/L) \n(L'/L) 

where a and 6 include derivatives like dY^(y, z)/dy\ y=0tZ=0 . With the intro- 
duction of the rescaling factor s = L'/L in equation (Bp, we can now evaluate 

Va- 

However an important point remains to be answered: how to estimate 
finite size dependence of v on lattice sizes Land// ? This can be achieved 
from large enough pairs of lattices L and V , with V > L, through the 
folowing expression for the partition function zeros: 

1 -rn(m^)/H$). do) 



vl,l> \ \ul{L) 

This equation defines our finite size estimators vl,l' from data in table 1. 
A second estimate can be obtained with the replacement \u\ — u c \ by its 
imaginary part Im(«5) in equation flTop. For large enough systems we have 
Re (ui) ~ u c , and both approaches should lead to the same result. 

In table 2 we present sequences of those two possible estimates i>l,sL as 
a function of the fixed rescaling factor s = 2, where the second column 
stands for the results of equation fllOD and the third one for the replacement 
\ul{sL)-u c \/\u\{L)-u c \ by \mu\{sL) / Im«?(L). 

As L increases, the values obtained by matching pairs of lattices are 
expected to converge to a limiting value. If we match our largest lattices 
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L = 24 and 32 we obtain from equation ( [3~Q|) v = 0.6260(66) whereas v = 
0.6292(70) is obtained by using only the imaginary part of the zeros. 

Looking at the values for the crossings (L, 2L) with (12, 24) and (16, 32) 
in the second column of table 2, we realize that our values for the real part of 
the corresponding zeros are not precise enough to obtain an ordered sequence 
towards its critical value v as L — > oo. On the other hand the estimate based 
on the imaginary part of zeros seems to preserve this trend. 

Since equations @ and fllCf ) are valid only for large enough lattice sizes, 
we have to discard the smallest values in table 2. For this reason we prefer 
not to evaluate equation @ by a multi-parameter fit. On the other hand, 
we note that this relation with V = sL, is an equation of type 

T(h) =T + ai h w + a 2 h 2w +•■■, (11) 

where we identify, 

V3 = ~W 

h = l/L 
T{h) = l/v L>2L 

(12) 

and T = \ jv the asymptotic limit for h — > 0. Hence, equation ( |TT1) is in the 
proper form to be analysed by the so called BST approximants, on which we 
elaborate in the next section. 



3.2 BST extrapolation 



Bulirsch and Stoer |H5| developed an algorithm to extrapolate a sequence h N , 
(N — 0, 1, 2, ...) converging to zero as iV — > oo. See also ref. [^] for a recent 
discussion on BST algorithm. 

The BST algorithm approximates tabulated data T(h) by a sequence of 
rational functions [pi], |2(| . The limiting value T is computed from a table of 
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recurrent relations defined from: 

T (AQ = Q 

T (7V) = T{h N ) (13) 

and 



T (N) = T (N+1) (T (N+1) _ T (N) 
m to— 1 ' \ m— 1 - L m— \, 



h \ - / rp(N+l) rp(N) 

n N \ I j _ TO — 1 —i TO-l I _ J 



(14) 



7 V -<■ TO— 1 1 TO — 2 

Here u> plays the role of a free parameter. If one defines e$ = 2(T^ +1 - ) — T$), 
it is expected — T\ < e$ in the limit % — >• oo. The above remark gives 
a criterion for choosing w as the value to minimize eW, in order to have 
a fast and reliable convergence. 

Our aim is to extrapolate the finite size sequence in table 2 by the BST 
algorithm. However, before proceeding with our analysis, we would like to 
explore first how the BST extrapolant approach performs for exact data. For 
this end, we come to the 2D Ising model. In ref. [^TJ the exact values for 



vl,2L are presented. We show in figures 2, 3, and 4 the BST extrapolants 
obtained from different sets of sequences T{hjsi). 

In figure 2 we show the BST estimates of the critical exponent z/ BST 
from a sequence for lattices of lengths L = h^ 1 , with h = 1/4, 1/6, 1/8, 
1/10, 1/12, 1/16, 1/20, 1/24 and 1/32. This figure presents a pole behavior 
at w ~ 1.580 for the extrapolated results (full lines) to h — > (L —>■ oo) for 
this sequence. In the neighbourhood of that pole we note the corresponding 
large values for the systematic error (dashed lines) according to the scale at 
the right hand side of that figure. Here we define the error as the difference 
between the extrapolated value T and a value of last but one interaction: 

rp rp{2) 

1 J m-1- 

We also observe that the known value v = 1 is obtained for w ~ 1 with 
error ~ 0. Moreover, it is remarkable that for this sequence of lengths h, z/ B 



8 



is weakly dependent on w. For instance, we obtain v with 0.1% precision, 
ue [0.999, 1.001] for a large range of w (we [0.1888, 1.5482]) before the pole. 

In figure 3 we restrict the available sequence to higher lattice sizes, h = 
1/12, 1/16, 1/20, 1/24 and 1/32. In that case the extrapolation with 0.1% 
precision is also compatible with a still large range of values for w: we 
[0.2816, 1.1969], before the pole at w ~ 1.2078. In figure 4, we restrict 
further the sequence to h = 1/16, 1/20, 1/24 and 1/32, and we obtain 
we [0.411, 1.3220], and a pole at w ~ 1.3871. 

Therefore, as we restrict our sequences to higher lattice sizes, the effects of 
the correction to scaling term become more pronounced, leading to a smaller 
range for acceptable values of w. This effect has a stronger counterpart in the 
criterium of minimum error: the acceptable range for w is actually narrower 
than stated above, mainly for figures 3 and 4. 

We are finally now at a point where we can use the BST algorithm to 
analyse our 3D Ising model data of table 2. Since we have to discard smaller 
lattices in order to accomplish with the expected validity conditions to derive 
equations (|) and (|H]), we are left with the sequence in the third column 
in table 2. Figure 5 presents our so obtained BST estimates for the critical 
exponent v from a sequence for lengths h = 1/6, 1/8, 1/12, and 1/16. Since 
our sequences for l/u L 2L are restricted to lattice sizes from L = 6 up to 16 
we take the four parameter fit estimate v = 0.6280 and its statistical error 
0.0015 as our input condition to find out w, as exemplified for the 2D Ising 
model. This statistical error leads to a range in w given by w = 0.745(74), 
marking the region of minimum systematic error. 

Finally, we remark that this approach yields a value for w which is in full 
agreement with recent results from various other and different approaches: 
from scaling relations of observables related to the magnetization || it was es- 
timated w = 0.87(9), from perturbative expansion at fixed D = 3 dimension 
followed w = 0.799(11), and w = 0.814(18) was estimated from e-expansion 
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4 Conclusions 



In summary, we have described a new calculation of the correction to 
scaling exponent w based on Monte Carlo multicanonical simulations and 
finite size scaling theory for the partition function zeros. A four-parameter 
fitting was done in order to find out the correlation length critical exponent 
v. Next, the result for v including the statistical error was used to obtain 
the acceptable values for the renormalization exponent y% = —w by means 
of the BST algorithm. This algorithm helped us to overcome the difficulties 
in the straight application of the multiparameter fit (|9|) to few data points 
and rather large lattice sizes. 

Our results for v and w are in good agreement with recently obtained esti- 
mates by Ballesteros et al. M as well as perturbative expansion calculations 
by Guida and Zinn- Justin J7|. 

It is tempting to assume that an accurate value for w could be pursued 
by increasing the significant precision of the complex partition function zeros 
of the 3D Ising model. This would account for a more precise calculation 
of vl,l>, by evaluating crossings between lattice sizes L and V . However as 
we have seen from figures 2 to 4, where we used exact values for vl,L', high 
precision in v does not necessarily lead to a smaller range in w. Hence, as 
an overall conclusion, we note that the large range of w is not just a matter 
of a lack of statistical precision but demonstrates that it is necessary to go 
to much larger lattice sizes. In particular, for 2D Ising model even L = 64 
seems to be not large enough. 
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Table Captions: 



Table 1: First partition function zeros for the 3D Ising model on a cubic 
lattice. 

Table 2: Sequence of estimates for the critical exponent vl,2L from pairs 
of lattices (L,2L). The third column is obtained by replacing |w5(2L) — 
u c \ I \ui(L) — u c \ by Im«5(2L) / ImtiJ(I) in equation fllPl) . 
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Table 1: First partition function zeros for the 3D Ising model on a cubic 
lattice. 



L 


NRUN 


NSWEEP 


Re(u?) 


Im(u?) 


Re (/??) 


Im (#) 


6 


2048 


100 000 


0.397586(18) 


0.045 435(17) 


0.228 964(11) 


0.028 445(11) 


8 


1024 


600 000 


0.402 723(11) 


0.028 596(10) 


0.226 748(07) 


0.017 722(06) 


12 


512 


900 000 


0.407018(12) 


0.014 925(12) 


0.224 557(07) 


0.009163(08) 


16 


512 


1200000 


0.408 814(11) 


0.009 422(11) 


0.223 557(07) 


0.005 761(07) 


24 


256 


1800 000 


0.410 341(14) 


0.004 935(12) 


0.222 674(09) 


0.003 006(07) 


32 


256 


1 800 000 


0.410 991(15) 


0.003 124(14) 


0.222 289(09) 


0.001900(08) 
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Table 2: Sequence of estimates for the critical exponent Vl,2L from pairs 
of lattices (L,2L). The third column is obtained by replacing \u\{2L) — 
u c \ I \u\(L) — u c \ by Imu5(2L) / Imw?(L) inequation (10). 



L 




VL,2L 


3 


0.60713(10) 


0.60916(11) 


4 


0.61986(13) 


0.61831(14) 


5 


0.62451(31) 


0.62181(33) 


6 


0.62577(44) 


0.62272(46) 


8 


0.62719(64) 


0.62429(67) 


12 


0.6278(14) 


0.6263(14) 


16 


0.6270(25) 


0.6279(26) 
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Figure Captions: 



Figure 1: Four parameter fit for \u^(L) — u c \ in function of L x l v in the 
range L = 3 — 32. The least square method gives v = 0.62853(35) with 
Q = 0.89. 

Figure 2: BST estimates of the critical exponent v (full lines) and sys- 
tematic error (dashed lines), as a function of the free parameter w, for 2D 
Ising model. The extrapolation is obtained from the finite size sequence with 
lattice sizes from L = 4 up to L = 32. The right side scale refers to the 
systematic error. 

Figure 3: As in figure 2, BST estimates of the critical exponent v (full lines) 
and systematic error (dashed lines) for 2D Ising model, with lattice sizes from 
L = 12 up to 32. The right side scale refers to the systematic error. 

Figure 4: As in figure 2, BST estimates of the critical exponent v (full lines) 
and systematic error (dashed lines) for 2D Ising model, with lattice sizes from 
L = 16 up to 32. The right side scale refers to the systematic error. 

Figure 5: BST estimates of the critical exponent v (full lines) and system- 
atic error (dashed lines), as a function of the free parameter w for 3D Ising 
model. The sequence is obtained with lattice sizes L — 6,8, 12 and 16. The 
right side scale refers to the systematic error. 
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Figure 1: Four parameter fit for \ui(L) — u c \ in function of L l l v in the 
range L = 3 — 32. The least square method gives v = 0.62853(35) with 
Q = 0.89. 
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Figure 2: BST estimates of the critical exponent v (full lines) and sys- 
tematic error (dashed lines), as a function of the free parameter w, for 2D 
Ising model. The extrapolation is obtained from the finite size sequence 
with lattice sizes from L = 4 up to L = 32. The right side scale refers to 
the systematic error. 



19 



CO 

> 



1.004 



1.002 - 



1.000 



0.998 



0.996 




0.002 



- 0.001 



0.000 



-0.001 



0.002 



CO 



Figure 3: As in figure 2, BST estimates of the critical exponent v (full 
lines) and systematic error (dashed lines) for 2D Ising model, with lattice 
sizes from L = 12 up to 32. The right side scale refers to the systematic 
error. 
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Figure 4: As in figure 2, BST estimates of the critical exponent v (full 
lines) and systematic error (dashed lines) for 2D Ising model, with lattice 
sizes from L = 16 up to 32. The right side scale refers to the systematic 
error. 
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Figure 5: BST estimates of the critical exponent v (full lines) and sys- 
tematic error (dashed lines), as a function of the free parameter w for 3D 
Ising model. The sequence is obtained with lattice sizes L — 6, 8, 12 and 
16. The right side scale refers to the systematic error. 
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